#=======================================================================================#
#
#  R code to replicate results reported in the paper titled
#  "A Forgotten Tie between Democracies and Nondemocracies in Asymmetric Cooperation"
#
#  Yasuki Kudo
#
#  Last modified: 2025-2-21
#=======================================================================================#

rm(list=ls()) # delete any objects in the memory

set.seed(100)

# packages
#=======================================================================================#
library(clusterSEs)
library(countrycode)
library(ggpubr)
library(haven)
library(miceadds)
library(marginaleffects)
library(modelsummary)
library(survival)
library(tidyverse)
library(xtable)
#devtools::install_github('IQSS/Zelig')
library(Zelig)
#install.packages("remotes")
#remotes::install_github("iqss/clarify")
library(clarify)
library(pollster)

# import data
#=======================================================================================#
orig.df <- read.csv("data/data_final.csv")

#=======================================================================================#
#=======================================================================================#
# Model 1 (Table 1 and Figures 1 and 2)
#=======================================================================================#
#=======================================================================================#

for(i in 1:9){ # seven models
  
  df <- orig.df
  
  if(i == 8){
    
    df$DemDem <- df$DemDem2
    df$DemNDem <- df$DemNDem2
    df$NDemDem <- df$NDemDem2
    df$NDemNDem <- df$NDemNDem2
    
  } 
  
  if(i == 9){
    
    df <- subset(df, mem1 != 2 & mem2 != 2)
    
  } 
  
  #### formulas ####
  
  # Main models (Table 1)
  
  if(i %in% c(1, 8, 9)){ 
    form <-  
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  # Table A2 and Figure A5
  
  if(i == 2){ # cinc 5% 
    form <- 
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither5) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  if(i == 3){ # cinc 20% 
    form <- 
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither20) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  if(i == 4){ # threats 20% 
    form <- 
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither20lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  if(i == 5){ # threats 40% 
    form <- 
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither40lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  if(i == 6){ # change in polity scores from Leeds & Savun (2007) 
    form <- 
      violate ~ 
      as.factor(DemDem) + as.factor(DemNDem) + as.factor(NDemDem) + as.factor(NDemNDem) + 
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regchangeeither) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  if(i == 7){ # main model (w/o as.factor()) for relogit
    form <-  
      violate ~ 
      DemDem + DemNDem + NDemDem + NDemNDem + 
      contig + capchangeeither + regimechange + chTEeither30lt + alform + colony +
      nomicoop + treaty + milinst + time + I(time^2) + I(time^3)
  }
  
  #### estimation ####
  if(!(i %in% c(7, 9))){
    
    mod <- glm(form, family=binomial(link='logit'), data = df) # standard SEs
    mod.cl <- miceadds::glm.cluster(form, family = binomial(link='logit'), data = df, cluster = "atopidphase") # cluster robust SEs
    vcov.mod.cl <- vcov(mod.cl)
    
  }
  
  if(i == 7){ # if relogit
    
    mod.zelig <- zelig(form, model = "relogit", data = df)
    mod <- from_zelig_model(mod.zelig)
    vcov.mod.cl <- vcovCL(mod, cluster = df$atopidphase)
    mod.cl <- coeftest(mod, vcov = vcov.mod.cl)
    
  }
  
  if(i == 9){
    
    mod <- glm(form, family=binomial(link='logit'), data = df, subset = (mem1 != 2 & mem2 != 2)) # standard SEs
    mod.cl <- miceadds::glm.cluster(form, family = binomial(link='logit'), data = subset(df, mem1 != 2 & mem2 != 2), cluster = "atopidphase") # cluster robust SEs
    vcov.mod.cl <- vcov(mod.cl)
    
  }
  
  # save vcov
  if(i == 1)vcov.mod.cl ->> vcov.mod.cl.main.all 
  if(i == 2)vcov.mod.cl ->> vcov.mod.cl.cinc5.all
  if(i == 3)vcov.mod.cl ->> vcov.mod.cl.cinc20.all
  if(i == 4)vcov.mod.cl ->> vcov.mod.cl.threats20.all
  if(i == 5)vcov.mod.cl ->> vcov.mod.cl.threats40.all
  if(i == 6)vcov.mod.cl ->> vcov.mod.cl.regime.all
  if(i == 7)vcov.mod.cl ->> vcov.mod.cl.relogit.all
  if(i == 8)vcov.mod.cl ->> vcov.mod.cl.pers.all
  if(i == 9)vcov.mod.cl ->> vcov.mod.cl.nonus.all
  
  # save models
  if(i == 1)mod ->> mod.main.all
  if(i == 2)mod ->> mod.cinc5.all
  if(i == 3)mod ->> mod.cinc20.all
  if(i == 4)mod ->> mod.threats20.all
  if(i == 5)mod ->> mod.threats40.all
  if(i == 6)mod ->> mod.regime.all
  if(i == 7)mod ->> mod.relogit.all
  if(i == 8)mod ->> mod.pers.all
  if(i == 9)mod ->> mod.nonus.all
  
  #### substantive effects ####
  if(i == 1){
    
    cat("\n")
    cat("\n")
    cat("######################## Substantive effects (Table 1) ########################")
    cat("\n")
    
    comp.mod <- 
      comparisons(mod, vcov = vcov.mod.cl, transform_pre = "ratio", conf_level = 0.95,
                  newdata = datagrid(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 0))
    
    cat("\n")
    cat("#### Dem-Dem ####")
    cat("\n")
    print((comp.mod[comp.mod$term == "DemDem", 4]-1)*100)
    cat("\n")
    cat("#### Dem-NDem ####")
    cat("\n")
    print((comp.mod[comp.mod$term == "DemNDem", 4]-1)*100)
    cat("\n")
    cat("#### NDem-Dem ####")
    cat("\n")
    print((comp.mod[comp.mod$term == "NDemDem", 4]-1)*100)
    cat("\n")
    cat("#### NDem-NDem ####")
    cat("\n")
    print((comp.mod[comp.mod$term == "NDemNDem", 4]-1)*100)
    cat("\n")
    cat("\n")
    
  }
  
  #### simulate coefficients ####
  if(i != 7)sim.mod <- sim(fit = mod, vcov = vcov.mod.cl)
  
  #### Figure for predicted prob. ####
  if(i == 1){
    
    time.max <- max(df$time, na.rm = T)
    time.min <- min(df$time, na.rm = T)
    
    pre.sym <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 0,
                        time = seq(time.min, time.max, 1)),
               verbose = FALSE)
    pre.sym <- t(pre.sym)
    pre.sym <- as.data.frame(pre.sym)
    p.mean.sym <- apply(pre.sym, 1, mean)
    pre.sym$p.mean <- p.mean.sym
    pre.sym$time <- seq(time.min, time.max, 1)
    pre.sym$type <- "sym"
    pre.sym <- select(pre.sym, p.mean, time, type)
    
    pre.demdem <- 
      sim_setx(sim.mod,
               x = list(DemDem = 1, DemNDem = 0, NDemDem = 0, NDemNDem = 0,
                        time = seq(time.min, time.max, 1)),
               verbose = FALSE)
    pre.demdem <- t(pre.demdem)
    pre.demdem <- as.data.frame(pre.demdem)
    p.mean.demdem <- apply(pre.demdem, 1, mean)
    pre.demdem$p.mean <- p.mean.demdem
    pre.demdem$time <- seq(time.min, time.max, 1)
    pre.demdem$type <- "demdem"
    pre.demdem <- select(pre.demdem, p.mean, time, type)
    
    pre.demndem <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0,
                        time = seq(time.min, time.max, 1)),
               verbose = FALSE)
    pre.demndem <- t(pre.demndem)
    pre.demndem <- as.data.frame(pre.demndem)
    p.mean.demndem <- apply(pre.demndem, 1, mean)
    pre.demndem$p.mean <- p.mean.demndem
    pre.demndem$time <- seq(time.min, time.max, 1)
    pre.demndem$type <- "demndem"
    pre.demndem <- select(pre.demndem, p.mean, time, type)
    
    pre.ndemdem <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 0, NDemDem = 1, NDemNDem = 0,
                        time = seq(time.min, time.max, 1)),
               verbose = FALSE)
    pre.ndemdem <- t(pre.ndemdem)
    pre.ndemdem <- as.data.frame(pre.ndemdem)
    p.mean.ndemdem <- apply(pre.ndemdem, 1, mean)
    pre.ndemdem$p.mean <- p.mean.ndemdem
    pre.ndemdem$time <- seq(time.min, time.max, 1)
    pre.ndemdem$type <- "ndemdem"
    pre.ndemdem <- select(pre.ndemdem, p.mean, time, type)
    
    pre.ndemndem <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 1,
                        time = seq(time.min, time.max, 1)),
               verbose = FALSE)
    pre.ndemndem <- t(pre.ndemndem)
    pre.ndemndem <- as.data.frame(pre.ndemndem)
    p.mean.ndemndem <- apply(pre.ndemndem, 1, mean)
    pre.ndemndem$p.mean <- p.mean.ndemndem
    pre.ndemndem$time <- seq(time.min, time.max, 1)
    pre.ndemndem$type <- "ndemndem"
    pre.ndemndem <- select(pre.ndemndem, p.mean, time, type)
    
    rbind(pre.sym, pre.demdem, pre.demndem, pre.ndemdem, pre.ndemndem) %>%
      ggplot(aes(time, p.mean, color = type)) +
      theme_bw() +
      geom_line(size = 1.5) +
      geom_hline(yintercept = 0, linewidth = 1, linetype = "dotted") +
      scale_color_manual(name="Alliance type", 
                         breaks  = c("sym", "ndemdem", "ndemndem", "demdem", "demndem"),
                         labels=c("Symm", "NDem-Dem", "Joint NDem", "Joint Dem", "Dem-NDem"),
                         values = c("1", "orange", "3", "4", "2")) +
      theme(text = element_text(size = 17),
            strip.text = element_text(size = 17),
            axis.text.y = element_text(colour = c('black')),
            axis.text.x = element_text(colour = c('black')),
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(2,20,2,20), 'mm'),
            legend.position = c(.8,.8),
            legend.direction="vertical",
            legend.key.width = unit(1,"cm"),
            legend.key.height = unit(.1,"cm"),
            legend.text = element_text(size = 15),
            legend.title = element_text(size = 15),) +
      ylab("P(Abrogation)") + xlab("Time since Formation (Years)") + ggtitle("Model 1") -> plot
    
    plot ->> pred.plot
    
    pre.demdem.2 <- 
      sim_setx(sim.mod,
               x = list(DemDem = 1, DemNDem = 0, NDemDem = 0, NDemNDem = 0), 
               verbose = FALSE)
    
    pre.demndem.2 <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0),
               verbose = FALSE)
    
    pre.ndemndem.2 <- 
      sim_setx(sim.mod,
               x = list(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 1),
               verbose = FALSE)
    
    cat("P(y|DemNDem) is lower than P(y|DemDem) by ")
    cat(((1-summary(pre.demndem.2)[1]/summary(pre.demdem.2)[1]))*100)
    cat("%")
    cat("\n")
    cat("\n")
    cat("P(y|DemNDem) is lower than P(y|NDemNDem) by ")
    cat(((1-summary(pre.demndem.2)[1]/summary(pre.ndemndem.2)[1]))*100)
    cat("%")
    cat("\n")
    cat("\n")
    cat("##########################################################################################")
  }
  
  #### Diff in predicted prob. ####
  if(i != 7){
    
    # DemNDem and Symm
    diff.sym <- sim_setx(sim.mod,
                         x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0),
                         x1 = list(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 0), 
                         verbose = FALSE)
    
    # DemNDem and DemDem
    diff.demdem <- sim_setx(sim.mod, 
                            x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0),
                            x1 = list(DemDem = 1, DemNDem = 0, NDemDem = 0, NDemNDem = 0), 
                            verbose = FALSE)
    
    # DemNDem and NDemDem
    diff.ndemdem <- sim_setx(sim.mod, 
                             x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0),
                             x1 = list(DemDem = 0, DemNDem = 0, NDemDem = 1, NDemNDem = 0), 
                             verbose = FALSE)
    
    # DemNDem and NDemNDem
    diff.ndemndem <- sim_setx(sim.mod, 
                              x = list(DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0),
                              x1 = list(DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 1), 
                              verbose = FALSE)
    
    p.diff <- 
      data.frame(fd = c(0, summary(diff.sym, level = .95)[3,1], summary(diff.demdem, level = .95)[3,1],
                        summary(diff.ndemdem, level = .95)[3,1], summary(diff.ndemndem, level = .95)[3,1]),
                 cil95 = c(0, summary(diff.sym, level = .95)[3,2], summary(diff.demdem, level = .95)[3,2],
                           summary(diff.ndemdem, level = .95)[3,2], summary(diff.ndemndem, level = .95)[3,2]),
                 cih95 = c(0, summary(diff.sym, level = .95)[3,3], summary(diff.demdem, level = .95)[3,3],
                           summary(diff.ndemdem, level = .95)[3,3], summary(diff.ndemndem, level = .95)[3,3]),
                 cil90 = c(0, summary(diff.sym, level = .9)[3,2], summary(diff.demdem, level = .9)[3,2],
                           summary(diff.ndemdem, level = .9)[3,2], summary(diff.ndemndem, level = .9)[3,2]),
                 cih90 = c(0, summary(diff.sym, level = .9)[3,3], summary(diff.demdem, level = .9)[3,3],
                           summary(diff.ndemdem, level = .9)[3,3], summary(diff.ndemndem, level = .9)[3,3]),
                 type = c("demndem", "sym", "demdem", "ndemdem", "ndemndem"))
    
    
  }
  
  if(i == 7){
    
    # DemNDem-Sym
    x <- setx(mod.zelig, DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0)
    x1 <- setx1(mod.zelig, DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 0)
    diff.sym <-  Zelig::sim(mod.zelig, x = x, x1 = x1, vcov = vcov.mod.cl)
    diff.sym <- as.list(diff.sym)
    diff.sym <- as.data.frame(diff.sym$sim.out$x1$fd[1])[,1]
    diff.sym.ci95 <- quantile(diff.sym, c(.025, .975))
    diff.sym.ci90 <- quantile(diff.sym, c(.05, .95))
    
    # DemNDem-DemDem
    x <- setx(mod.zelig, DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0)
    x1 <- setx1(mod.zelig, DemDem = 1, DemNDem = 0, NDemDem = 0, NDemNDem = 0)
    diff.demdem <-  Zelig::sim(mod.zelig, x = x, x1 = x1, vcov = vcov.mod.cl)
    diff.demdem <- as.list(diff.demdem)
    diff.demdem <- as.data.frame(diff.demdem$sim.out$x1$fd[1])[,1]
    diff.demdem.ci95 <- quantile(diff.demdem, c(.025, .975))
    diff.demdem.ci90 <- quantile(diff.demdem, c(.05, .95))
    
    # DemNDem-NDemDem
    x <- setx(mod.zelig, DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0)
    x1 <- setx1(mod.zelig, DemDem = 0, DemNDem = 0, NDemDem = 1, NDemNDem = 0)
    diff.ndemdem <-  Zelig::sim(mod.zelig, x = x, x1 = x1, vcov = vcov.mod.cl)
    diff.ndemdem <- as.list(diff.ndemdem)
    diff.ndemdem <- as.data.frame(diff.ndemdem$sim.out$x1$fd[1])[,1]
    diff.ndemdem.ci95 <- quantile(diff.ndemdem, c(.025, .975))
    diff.ndemdem.ci90 <- quantile(diff.ndemdem, c(.05, .95))
    
    # DemNDem-NDemNDem
    x <- setx(mod.zelig, DemDem = 0, DemNDem = 1, NDemDem = 0, NDemNDem = 0)
    x1 <- setx1(mod.zelig, DemDem = 0, DemNDem = 0, NDemDem = 0, NDemNDem = 1)
    diff.ndemndem <-  Zelig::sim(mod.zelig, x = x, x1 = x1, vcov = vcov.mod.cl)
    diff.ndemndem <- as.list(diff.ndemndem)
    diff.ndemndem <- as.data.frame(diff.ndemndem$sim.out$x1$fd[1])[,1]
    diff.ndemndem.ci95 <- quantile(diff.ndemndem, c(.025, .975))
    diff.ndemndem.ci90 <- quantile(diff.ndemndem, c(.05, .95))
    
    
    p.diff <- 
      data.frame(fd = c(0, mean(diff.sym.ci95), mean(diff.demdem.ci95), mean(diff.ndemdem.ci95), mean(diff.ndemndem.ci95)),
                 cil95 = c(0, diff.sym.ci95[1], diff.demdem.ci95[1], diff.ndemdem.ci95[1], diff.ndemndem.ci95[1]),
                 cih95 = c(0, diff.sym.ci95[2], diff.demdem.ci95[2], diff.ndemdem.ci95[2], diff.ndemndem.ci95[2]),
                 cil90 = c(0, diff.sym.ci90[1], diff.demdem.ci90[1], diff.ndemdem.ci90[1], diff.ndemndem.ci90[1]),
                 cih90 = c(0, diff.sym.ci90[2], diff.demdem.ci90[2], diff.ndemdem.ci90[2], diff.ndemndem.ci90[2]),
                 type = c("demndem", "sym", "demdem", "ndemdem", "ndemndem"))
    
  }
  
  if(i == 1)title <- ggtitle("Model 1")
  if(i == 2)title <- ggtitle("CINC 5%")
  if(i == 3)title <- ggtitle("CINC 20%")
  if(i == 4)title <- ggtitle("Threats 20%")
  if(i == 5)title <- ggtitle("Threats 40%")
  if(i == 6)title <- ggtitle("L&S's (07) regime change variable")
  if(i == 7)title <- ggtitle("Rare event logit")
  if(i == 8)title <- ggtitle("Personalistic authority")
  if(i == 9)title <- ggtitle("Non-US alliances")
  
  p.diff %>%
    ggplot(aes(NULL)) +
    theme_bw() +
    geom_linerange(aes(x = type, y = fd, ymax = cih90, ymin = cil90), linewidth = 3.5, color = "gray50") +
    geom_pointrange(aes(x = type, y = fd, ymax = cih95, ymin = cil95), size = 1, linewidth = 1.5) +
    geom_text(aes(x = type, y = fd, label = round(fd, digits = 3), vjust = -0.8), size = 5.5) +
    geom_hline(yintercept = 0, linetype = "dotted", size = .5) +
    scale_x_discrete(limit = c("ndemndem", "ndemdem", "demdem", "sym", "demndem"),
                     label = c("Joint NDem", "NDem-Dem", "Joint Dem", "Symm", "Dem-NDem \n (Baseline)")) +
    xlab("") + ylab("Diff. in P(Abrogation)") + title  +
    scale_color_identity() +
    coord_flip() +
    theme(text = element_text(size = 17),
          axis.text.y = element_text(size = 17, colour = c('black')),
          axis.text.x = element_text(colour = c('black')),
          strip.text = element_text(size = 17),
          plot.margin=unit(c(2,2,2,2), 'mm')) -> plot
  
  if(i == 1)plot ->> fd.plot.main
  if(i == 2)plot ->> fd.plot.cinc5
  if(i == 3)plot ->> fd.plot.cinc20
  if(i == 4)plot ->> fd.plot.threats20
  if(i == 5)plot ->> fd.plot.threats40
  if(i == 6)plot ->> fd.plot.regime
  if(i == 7)plot ->> fd.plot.relogit
  if(i == 8)plot ->> fd.plot.pers
  if(i == 9)plot ->> fd.plot.nonus
  
}

#### tables ####
# table 1
coef_map <- c(`as.factor(DemDem)1` = "Joint Dem Asymm",
              `as.factor(DemNDem)1` = "Dem-Nondem Asymm",
              `as.factor(NDemDem)1` = "Nondemo-Dem Asymm",
              `as.factor(NDemNDem)1` = "Joint Nondem Asymm",
              `as.factor(majdem)1` = "Asymm Alliance w/ Democratic Major Power",
              `as.factor(majndem)1` = "Asymm Alliance w/ Nondemocratic Major Power",
              `as.factor(mindem)1` = "Asymm Alliance w/ Democratic Minor Power",
              `as.factor(minndem)1` = "Asymm Alliance w/ Nondemocratic Minor Power",
              `as.factor(contig)1` = "Contiguity",
              `as.factor(capchangeeither)1` = "Change in National Capabilities",
              `as.factor(capchangeeither5)1` = "Change in National Capabilities",
              `as.factor(capchangeeither20)1` = "Change in National Capabilities",
              `as.factor(regchangeeither)1` = "Change in Polity",
              `as.factor(regimechange)1` = "Recent Regime Change",
              `as.factor(chTEeither20lt)1` = "Change in External Threats",
              `as.factor(chTEeither30lt)1` = "Change in External Threats",
              `as.factor(chTEeither40lt)1` = "Change in External Threats",
              `as.factor(alform)1` = "Alliance Formation",
              `as.factor(nomicoop)1` = "Non-Miilitary Cooperation",
              `as.factor(treaty)1` = "Treaty Ratification",
              `DemDem` = "Joint Dem Asymm",
              `DemNDem` = "Dem-Nondem Asymm",
              `NDemDem` = "Nondemo-Dem Asymm",
              `NDemNDem` = "Joint Nondem Asymm",
              `majdem` = "Asymm Alliance w/ Democratic Major Power",
              `majndem` = "Asymm Alliance w/ Nondemocratic Major Power",
              `mindem` = "Asymm Alliance w/ Democratic Minor Power",
              `minndem` = "Asymm Alliance w/ Nondemocratic Minor Power",
              `ndem` = "Number of Democracies",
              `contig` = "Contiguity",
              `capchangeeither` = "Change in National Capabilities",
              `capchangeeither5` = "Change in National Capabilities",
              `capchangeeither20` = "Change in National Capabilities",
              `regchangeeither` = "Change in Polity",
              `regimechange` = "Recent Regime Change",
              `chTEeither20lt` = "Change in External Threats",
              `chTEeither30lt` = "Change in External Threats",
              `chTEeither40lt` = "Change in External Threats",
              `alform` = "Alliance Formation",
              `nomicoop` = "Non-Miilitary Cooperation",
              `treaty` = "Treaty Ratification",
              `milinst` = 'Military Coordination',
              `as.factor(colony)1` = 'Colony',
              `colony` = 'Colony')

# Table 1
modelsummary(list("Model 1" = mod.main.all), 
             vcov = list(vcov.mod.cl.main.all), 
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, align = "ld", gof_map = c("nobs"),
             title = "Logistic Regression of Alliance Abrogation")
logLik(mod.main.all)

# Table A3
model.robust.all <- 
  list("CINC 5%" = mod.cinc5.all, "CINC 20%" = mod.cinc20.all, "Threats 20%" = mod.threats20.all, 
       "Threats 40%" = mod.threats40.all, "L&S07" = mod.regime.all, "Relogit" = mod.relogit.all,
       "Non-US alliances" = mod.nonus.all)

vcov.robust.all <- 
  list(vcov.mod.cl.cinc5.all, vcov.mod.cl.cinc20.all, vcov.mod.cl.threats20.all,
       vcov.mod.cl.threats40.all, vcov.mod.cl.regime.all, vcov.mod.cl.relogit.all,
       vcov.mod.cl.nonus.all)

modelsummary(model.robust.all, vcov = vcov.robust.all, 
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, align = "lddddddd", gof_map = c("nobs"),
             title = "Robustness checks", 
             #output = "latex"
             )
logLik(mod.cinc5.all)
logLik(mod.cinc20.all)
logLik(mod.threats20.all)
logLik(mod.threats40.all)
logLik(mod.regime.all)
logLik(mod.relogit.all)
logLik(mod.nonus.all)

# Table A6
modelsummary(mod.pers.all, vcov = vcov.mod.cl.pers.all, 
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, align = "ld", gof_map = c("nobs"),
             title = "Robustness checks", 
             #output = "latex"
             )
logLik(mod.pers.all)


#### figures ####
ggarrange(fd.plot.cinc5, fd.plot.cinc20,
          fd.plot.threats20, fd.plot.threats40,
          fd.plot.regime, fd.plot.relogit,
          fd.plot.nonus) -> fd.plot.robust

ggsave("SSQ/paper/figures/pred_plot.pdf", pred.plot, width = 8, height = 6)         # Figure 1
ggsave("SSQ/paper/figures/fd_plot_main.pdf", fd.plot.main, width = 8, height = 4)   # Figure 2
ggsave("SSQ/paper/figures/fd_robust.pdf", fd.plot.robust, width = 24, height = 12)  # Figure A5
ggsave("SSQ/paper/figures/fd_plot_pers.pdf", fd.plot.pers, width = 8, height = 4)   # Figure A6


#=======================================================================================#
#=======================================================================================#
# Monadic effects of a major power and a minor power
#=======================================================================================#
#=======================================================================================#

#### Model estimation and testing diff in coefficients
for(i in 1:4){
  
  df <- orig.df
  
  cat("\n")
  cat("\n")
  cat("\n")
  cat("##################################################################")
  cat("\n")
  if(i == 1)print("#### Model 2 (Table A4) ####")
  if(i == 2)print("#### Model 3 (Table A4) ####")
  if(i == 3)print("#### Model 2 (Table A5) ####")
  if(i == 4)print("#### Model 3 (Table A5)  ####")
  cat("##################################################################")
  cat("\n")
  
  #### formulas ####
  # Model 2
  if(i == 1){
    form <- 
      violate ~ 
      as.factor(majdem) + as.factor(majndem) +
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  # Model 3 
  if(i == 2){
    form <- 
      violate ~ 
      as.factor(mindem) + as.factor(minndem) +
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  # Model 2 (Table A3)
  if(i == 3){
    form <- 
      violate ~ 
      as.factor(majdem) + as.factor(majndem) +  ndem +
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  # Model 3 (Table A3)
  if(i == 4){
    form <- 
      violate ~ 
      as.factor(mindem) + as.factor(minndem) + ndem +
      as.factor(contig) + as.factor(capchangeeither) + as.factor(regimechange) + as.factor(chTEeither30lt) + as.factor(alform) + as.factor(colony) +
      as.factor(nomicoop) + as.factor(treaty) + milinst + time + I(time^2) + I(time^3)
  }
  
  
  #### estimation ####
  mod <- glm(form, family=binomial(link='logit'), data = df) # standard SEs
  mod.cl <- miceadds::glm.cluster(form, family=binomial(link='logit'), data = df, cluster = "atopidphase") # cluster robust SEs
  vcov.cl <- vcov(mod.cl)
  
  
  #### test diffs in coefficients ####
  
  cat("\n")
  cat("\n")
  cat("######################## Testing diff. in coefficients ########################")
  cat("\n")
  cat("\n")
  
  if(i %in% c(1, 3))print(car::linearHypothesis(mod, vcov = vcov.cl, c("as.factor(majdem)1 = as.factor(majndem)1")))
  if(i %in% c(2, 4))print(car::linearHypothesis(mod, vcov = vcov.cl, c("as.factor(mindem)1 = as.factor(minndem)1")))
  
  if(i == 1)mod ->> mod.maj
  if(i == 2)mod ->> mod.min
  if(i == 3)mod ->> mod.maj.ndem
  if(i == 4)mod ->> mod.min.ndem
  
  if(i == 1)vcov.cl ->> vcov.maj
  if(i == 2)vcov.cl ->> vcov.min
  if(i == 3)vcov.cl ->> vcov.maj.ndem
  if(i == 4)vcov.cl ->> vcov.min.ndem
  
  #### substantive effects ####
  cat("\n")
  cat("\n")
  cat("######################## Percent changes ########################")
  cat("\n")
  
  if(i %in% c(1, 3)){
    comp <- 
      comparisons(mod, vcov = vcov.cl, transform_pre = "ratio",
                  newdata = datagrid(majdem = 0, majndem = 0))
    
    cat("\n")
    cat("#### Asymm Alliance w/ Democratic Major Power ####")
    cat("\n")
    print((comp[comp$term == "majdem", 4]-1)*100)
    cat("\n")
    cat("#### Asymm Alliance w/ Nondemocratic Major Power ####")
    cat("\n")
    print((comp[comp$term == "majndem", 4]-1)*100)
    
  }
  
  if(i %in% c(2, 4)){
    comp <- 
      comparisons(mod, vcov = vcov.cl, transform_pre = "ratio",
                  newdata = datagrid(mindem = 0, minndem = 0))
    cat("\n")
    cat("#### Asymm Alliance w/ Democratic Minor Power ####")
    cat("\n")
    print((comp[comp$term == "mindem", 4]-1)*100)
    cat("\n")
    cat("#### Asymm Alliance w/ Nondemocratic Minor Power ####")
    cat("\n")
    print((comp[comp$term == "minndem", 4]-1)*100)
    
  }
  
}

# Table A4
model.monad <- list("Maj pow" = mod.maj, "Min pow" = mod.min)
vcov.monad <- list(vcov.maj, vcov.min)

modelsummary(models = model.monad, vcov = vcov.monad, 
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, align = "ldd", gof_map = c("nobs", "logLik"),
             title = "Monadic effects"
             #output = "latex",
             )

# Table A5
model.monad.ndem <- list("Maj pow" = mod.maj.ndem, "Min pow" = mod.min.ndem)
vcov.monad.ndem <- list(vcov.maj.ndem, vcov.min.ndem)

modelsummary(models = model.monad.ndem, vcov = vcov.monad.ndem, 
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, align = "ldd", gof_map = c("nobs", "logLik"),
             title = "Monadic effects",
             #output = "latex"
             )

#=======================================================================================#
#=======================================================================================#
# Cox proportional hazards model
#=======================================================================================#
#=======================================================================================#

### formulas ### 
# Model 1
formula.m1.surv <- 
  Surv(dur1, dur2, violate) ~
  DemDem + DemNDem + NDemDem + NDemNDem +
  contig + capchangeeither + regimechange + chTEeither30lt + 
  alform + colony + nomicoop + treaty + milinst

# Model 2
formula.m2.surv <- 
  Surv(dur1, dur2, violate) ~
  majdem + majndem +
  contig + capchangeeither + regimechange + chTEeither30lt + 
  alform + colony + nomicoop + treaty + milinst

# Model 3
formula.m3.surv <- 
  Surv(dur1, dur2, violate) ~
  mindem + minndem +
  contig + capchangeeither + regimechange + chTEeither30lt + 
  alform + colony + nomicoop + treaty + milinst

tidy_custom.coxph <- function(x) {
  s <- summary(x)$coefficients
  data.frame(
    term = row.names(s),
    robust.se = s[, "robust se", drop = FALSE])
}

### Model 1 ###
# estimation
summary(mod1.cox <- coxph(formula.m1.surv, data = orig.df, id = atopidphase, cluster = atopidphase))
vcov1.cox <- vcov(mod1.cox)
logLik(mod1.cox)

# substantive effects
summary(comparisons(mod1.cox, comparison = "ratio", variables = list(DemDem = c(0, 1))))
summary(comparisons(mod1.cox, comparison = "ratio", variables = list(DemNDem = c(0, 1))))
summary(comparisons(mod1.cox, comparison = "ratio", variables = list(NDemDem = c(0, 1))))
summary(comparisons(mod1.cox, comparison = "ratio", variables = list(NDemNDem = c(0, 1))))

### Model 2 ###
# estimation
summary(mod2.cox <- coxph(formula.m2.surv, data = orig.df, id = atopidphase, cluster = atopidphase))
vcov2.cox <- vcov(mod2.cox)
logLik(mod2.cox)

# Diff in coeffs
print(car::linearHypothesis(mod2.cox, vcov = vcov2.cox, c("majdem = majndem")))

# substantive effects
summary(comparisons(mod2.cox, comparison = "ratio", variables = list(majdem = c(0, 1))))
summary(comparisons(mod2.cox, comparison = "ratio", variables = list(majndem = c(0, 1))))


### Model 3 ###
# estimation
summary(mod3.cox <- coxph(formula.m3.surv, data = orig.df, id = atopidphase, cluster = atopidphase))
vcov3.cox <- vcov(mod3.cox)
logLik(mod3.cox)

# Diff in coeffs
print(car::linearHypothesis(mod3.cox, vcov = vcov3.cox, c("mindem = minndem")))

# substantive effects
summary(comparisons(mod3.cox, comparison = "ratio", variables = list(mindem = c(0, 1))))
summary(comparisons(mod3.cox, comparison = "ratio", variables = list(minndem = c(0, 1))))

#### tables ####
# Table A2
modelsummary(list("Model 1 (Cox model)" = mod1.cox), 
             statistic = "robust.se",
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, gof_map = "nobs",
             title = "Model 1 using Cox proportional hazards model")

# Models 2 and 3 with the COX model, not reported in the paper
modelsummary(list("Model 2 (Cox model)" = mod2.cox, "Model 3 (Cox model)" = mod3.cox), 
             statistic = "robust.se",
             stars = c('+' = .1, '*' = .05, '**' = .01),
             coef_map = coef_map, gof_map = "nobs",
             title = "Models 2 and 3 using Cox proportional hazards model")


#=======================================================================================#
#=======================================================================================#
# information shown in Appendix
#=======================================================================================#
#=======================================================================================#

# Figure A.1
time <- seq(1, 51, 1)

p_vio <- table(orig.df$time, orig.df$violate)[,2]/table(orig.df$time, orig.df$violate)[,1]

ggplot(NULL) +
  theme_bw() +
  geom_vline(xintercept = c(quantile(orig.df$time, probs = c(0.05, 0.5, 0.95))),
             size = 1, linetype = "dotted") +
  geom_hline(yintercept = 0, size = 1, colour = "gray") +
  geom_point(aes(time, p_vio), size = 3) +
  geom_smooth(aes(time, p_vio), se = F) +
  geom_point(size = 6) +
  xlab("Time since Formation (Years)") + 
  ylab("Proportion of Abrogation") +
  theme(text = element_text(size = 17),
        plot.title = element_text(size = 17),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        plot.margin=unit(c(3,3,3,3), 'mm')) -> p_timevio

p_timevio

# Figure A.2
ggplot(NULL) +
  theme_bw() +
  geom_histogram(aes(orig.df$time), binwidth = 1, color = "black", fill = "gray") +
  geom_vline(xintercept = c(quantile(orig.df$time, probs = c(0.05, 0.5, 0.95))),
             size = 1, linetype = "dotted") +
  xlab("Time since Formation (Years)") + 
  ylab("Frequency") +
  theme(text = element_text(size = 17),
        plot.title = element_text(size = 17),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        plot.margin=unit(c(3,3,3,3), 'mm')) -> hist_timevio

hist_timevio

# Figure A.3
ggplot() +
  theme_bw() +
  geom_bar(aes(as.factor(alliance)), subset(orig.df, is.na(alliance) == F), color="black", fill="gray",  na.rm = T) +
  theme(text = element_text(size = 17),
        strip.text = element_text(size = 17),
        strip.background = element_rect(colour = 'black',
                                        fill = 'white', 
                                        linetype='solid'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 15,
                                   color = "black"
        ),
        axis.text.y = element_text(size = 17,
                                   color = "black"),
        plot.margin=unit(c(3,3,3,3), 'mm')) +
  coord_flip() +
  scale_x_discrete(limit = factor(c(4, 3, 2, 1, 0)),
                   label = c( "Nondem.-Nondem. \n (70)",
                              "Nondem.-Dem. \n (14)",
                              "Dem.-NonDem. \n (48)", 
                              "Dem.-Dem. \n (22)", 
                              "Sym. Alliance \n (139)")) +
  xlab("Alliance Type") +  ylab("Frequency") -> fre_alliance

fre_alliance

# Figure A.4
orig.df_propvio <- 
  tibble(n = table(orig.df$alliance, orig.df$violate)[,1],
         n_vio = table(orig.df$alliance, orig.df$violate)[,2],
         p = n_vio/n,
         ciu = p+1.96*sqrt((p*(1-p))/n),
         cil = p-1.96*sqrt((p*(1-p))/n),
         type = c(0, 1, 2, 3, 4)) 


ggplot(NULL) +
  theme_bw() +
  geom_hline(yintercept = 0, size = 1, linetype = "dotted") +
  geom_linerange(aes(x = as.factor(type), ymax = ciu, ymin = cil), orig.df_propvio,
                 size = 1) +
  geom_point(aes(x = as.factor(type), y = p), orig.df_propvio, size = 3.5) +
  theme(text = element_text(size = 17),
        strip.text = element_text(size = 17),
        strip.background = element_rect(colour = 'black',
                                        fill = 'white', 
                                        linetype='solid'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 15,
                                   color = "black"),
        axis.text.y = element_text(size = 17,
                                   color = "black"),
        plot.margin=unit(c(3,3,3,3), 'mm')) +
  coord_flip() +
  scale_x_discrete(limit = factor(c(4, 3, 2, 1, 0)),
                   label = c( "Nondem.-Nondem.",
                              "Nondem.-Dem.",
                              "Dem.-NonDem.", 
                              "Dem.-Dem.", 
                              "Sym. Alliance")) +
  ylab("Proportion of Abrogation") + xlab("Alliance Type") -> plot_p_vio

# save figures
ggsave("SSQ/paper/figures/hist_timevioPlot.pdf", hist_timevio, width = 8, height = 7) # Figure A1
ggsave("SSQ/paper/figures/p_timevioPlot.pdf", p_timevio, width = 8, height = 7)      # Figure A2
ggsave("SSQ/paper/figures/fre_alliancePlot.pdf", fre_alliance, width = 10, height = 5)  # Figure A3
ggsave("SSQ/paper/figures/alliancevioPlot.pdf", plot_p_vio, width = 10, height = 5)     # Figure A4

## List of alliance Abrogation
orig.df %>%
  group_by(atopidphase, alliance) %>%
  summarise(n = n()) -> n.dyad

table(n.dyad$alliance)

orig.df %>%
  filter(violate == 1 & !is.na(alliance)) %>%
  mutate(Democracy1 = ifelse(dem1 == 1, "Yes", "No"),
         Democracy2 = ifelse(dem2 == 1, "Yes", "No"),
         majpow1 = ifelse(majpow1 == 1, "Yes", "No"),
         majpow2 = ifelse(majpow2 == 1, "Yes", "No"),) %>%
  left_join(data.frame(name1 = codelist_panel$country.name.en, year = codelist_panel$year, mem1 = codelist_panel$cown),
            by = c("mem1", "year")) %>%
  left_join(data.frame(name2 = codelist_panel$country.name.en, year = codelist_panel$year, mem2 = codelist_panel$cown),
            by = c("mem2", "year")) %>%
  arrange(year, dem1, dem2) %>%
  dplyr::select(atopidphase, year, name1, name2, Democracy1, Democracy2, majpow1, majpow2) -> list.vio

xtable(list.vio, digit = 0, caption = "List of Alliancen Abrogation", label = "list_vio") # Table A1




